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1 Introduction 

Plasmas, which are gases of charged particles, and charged particle beams 
can be described by a distribution function f(t,x,v) dependent on time t, 
on position x and on velocity v. The function / represents the probability of 
presence of a particle at position (x, v) in phase space at time t. It satisfies 
the so-called Vlasov equation 

%+v-V x f + F(t,x,v)-V v f = 0. (1) 
at 

The force field F(t, x, v) consists of applied and self-consistent electric and 
magnetic fields: 

F = —(E se if + E app + v x (B se if + B app )), 

where m represents the mass of a particle and q its charge. The self-consistent 
part of the force field is solution of Maxwell's equations 

1 <9E p 

--l^ + VxB = ^j' V-E=-S 

c z ot e 

?|VxE = 0, V-B = 0. 

ot 

The coupling with the Vlasov equation results from the source terms p 
and j such that: 

p(t,x)=q f(t,x,v)dv, i = q I f(t,x,v)vdv. 

We then obtain the nonlinear Vlasov-Maxwell equations. In some cases, when 
the field are slowly varying the magnetic field becomes negligible and the 
Maxwell equations can be replaced by the Poisson equation where: 

E se tf(t,x) = -V x <f>(t,x), -e Q A x (j) = p. (2) 



The numerical resolution of the Vlasov equation is usually performed 
by particle methods (PIC) which consist in approximating the plasma by a 
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finite number of particles. The trajectories of these particles are computed 
from the characteristic curves given by the Vlasov equation, whereas self- 
consistent fields are computed on a mesh of the physical space. This method 
allows to obtain satisfying results with a few number of particles. However, 
it is well known that, in some cases, the numerical noise inherent to the 
particle method becomes too important to have an accurate description of 
the distribution function in phase space. Moreover, the numerical noise only 
decreases in yN, when the number of particles N is increased. To remedy to 
this problem, methods discretizing the Vlasov equation on a mesh of phase 
space have been proposed. A review of the main methods for the resolution 
of the Vlasov equation is given in these proceedings [5] . 

The major drawback of methods using a uniform and fixed mesh is that 
their numerical cost is high, which makes them rather inefficient when the 
dimension of phase-space grows. For this reason we are investigating here 
a method using an adaptive mesh. The adaptive method is overlayed to a 
classical semi-Lagrangian method which is based on the conservation of the 
distribution function along characteristics. Indeed, this method uses two steps 
to update the value of the distribution function at a given mesh point. The 
first one consists in following the characteristic ending at this mesh point 
backward in time, and the second one in interpolating its value there from 
the old values at the surrounding mesh points. Using the conservation of the 
distribution function along the characteristics this will yield its new value 
at the given mesh point. This idea was originally introduced by Cheng and 
Knorr [2j along with a time splitting technique enabling to compute exactly 
the origin of the characteristics at each fractional step. In the original method, 
the interpolation was performed using cubic splines. This method has since 
been used extensively by plasma physicists (see for example [H [6] and the ref- 
erences therein). It has then been generalized to the frame of semi-Lagrangian 
methods by E. Sonnendriicker et al. [8]. This method has also been used to 
investigate problems linked to the propagation of strongly nonlinear heavy 
ion beams [9]- 

In the present work, we have chosen to introduce a phase-space mesh 
which can be refined or derefined adaptively in time. For this purpose, we 
use a technique based on multiresolution analysis which is in the same spirit 
as the methods developed in particular by S. Bertoluzza pQ, A. Cohen et al. 
[3] and M. Griebel and F. Koster [7]. We represent the distribution function 
on a wavelet basis at different scales. We can then compress it by eliminat- 
ing coefficients which are small and accordingly remove the associated mesh 
points. Another specific feature of our method is that we use an advection in 
physical and velocity space forward in time to predict the useful grid points 
for the next time step, rather than restrict ourselves to the neighboring points. 
This enables us to use a much larger time step, as in the semi-Lagrangian 
method the time step is not limited by a Courant condition. Once the new 
mesh is predicted, the semi-Lagrangian methodology is used to compute the 
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new values of the distribution function at the predicted mesh points, using 
an interpolation based on the wavelet decomposition of the old distribution 
function. The mesh is then refined again by performing a wavelet transform, 
and eliminating the points associated to small coefficients. 

This paper is organized as follows. In section [2l we recall the tools of 
multiresolution analysis which will be needed for our method, precizing what 
kind of wavelets seem to be the most appropriate in our case. Then, we 
describe in section [3] the algorithm used in our method, first for the non 
adaptive mesh case and then for the adaptive mesh case. Finally we present 
a few preliminary numerical results. 

2 Multiresolution analysis 

The semi-Lagrangian method consists mainly of two steps, an advection step 
and an interpolation step. The interpolation part is performed using for ex- 
ample a Lagrange interpolating polynomial on a uniform grid. Thus interpo- 
lating wavelets provide a natural way to extend this procedure to an adaptive 
grid in the way we shall now shortly describe. 

For simplicity, we shall restrict our description to the ID case of the whole 
real line. It is straightforward to extend it to periodic boundary conditions 
and it can also be extended to an interval with Dirichlct boundary conditions. 
The extension to higher dimension is performed using a tensor product of 
wavelets and will be addressed at the end of the section. 

For any value of j 6 Z, we consider a uniform grid G J of step 2" J . The 
grid points are located at x J k = k2~ 3 . This defines an infinite sequence of 
grids that we denote by (Gj)j^z, and j will be called the level of the grid. 

In order to go from one level to the next or the previous, we define a pro- 
jection operator and a prediction operator. Consider two grid levels Gj and 
Gj+i and discrete values (of a function) denoted by {d k )k^z and {c? k +1 )kez. 
Even though we use the same index k for the grid points in the two cases, 
there are of course twice as many points in any given interval on Gj+i as on 
Gj. Using the terminology in [3], we then define the projection operator 

Pj+i '■ G\H-i — * ( ' r 
which is merely a restriction operator, as well as the prediction operator 

such that e|£ = c£, 

4fc+l = P2N+l(%lk+l)> 
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where P2N+1 stands for the Lagrange interpolation polynomial of odd degree 
2N + 1 centered at the point (a^fe+i)- 

Using the just defined prediction operator, we can construct on Gj a 
subspace of L 2 (R) that we shall denote by Vj, a basis of which being given 
by {<p k )kei, such that <p k (x k ,) = 8 kk i where S kk i is the Kronecker symbol. The 
value of ipL at any point of the real line is then obtained by applying, possibly 
an infinite number of times, the prediction operator. 

In the wavelets terminology the p? k are called scaling functions. We shall 
also denote by p — <ps®. Let us notice that 

<fi(x)=<p(2>x-k). 

It can be easily verified that the scaling functions satisfy the following prop- 
erties: 

— Compact support: the support of ip is included in [— 2N — 1, 2N + 1]. 

— Interpolation: by construction (p(x) is interpolating in the sense that <^(0) = 
1 and ip(k) = if k ^ 0. 

— Polynomial representation: all polynomials of degree less or equal to 2N+ 1 
can be expressed exactly as linear combinations of the p k . 

— Change of scale: the p at a given scale can be expressed as a linear combi- 
nation of the ip at the scale immediately below: 

2N+1 

p(x) = ^ hip(2x — I). 

-2N-1 

Moreover the sequence of spaces {Vj )j e z defines a multircsolution analysis 
of L 2 (R), i.e. it satisfies the following properties: 

— . . . C VLi C Vp C V1C...CKC...C L 2 (R). 

— r\Vj = {0}, UVj = L 2 (R). 

— feV j ^f(2-)V j+1 . 

— 3p (scaling function) such that {p{x — fc)}fcez is a basis of Vo and {<p? k = 
2^l' 2 ip{2? x — fc)}fcgz is a basis of Vj. 

As Vj C Vj+i, there exists a supplementary of Vj in Vj + i that we shall 
call the detail space and denote by Wj : 

Vj +l =Vj@Wj. 

The construction of Wj can be made in the following way: an element of 
Vj+\ is characterized by the sequence(c^ +1 )fc e z and by construction we have 
c\ = 4ft 1 - Thus, if we define d\ = 4tfi _ ^JV+iO^fc+i): where P 2 n+i is 
the Lagrange interpolation polynomial by which the value of an element of 
Vj at the point (ar^fe+i) can be computed, d? k represents exactly the difference 
between the value in Vj + i and the value predicted in Vj. Finally, any element 
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of Vj + \ can be characterized by the two sequences {cjjk of values in Vj 
and (d J k )k of details in Wj. Moreover this strategy for constructing Wj is 
particularly interesting for adaptive refinement as d? k will be small at places 
where the prediction from Vj is good and large elsewhere, which gives us a 
natural refinement criterion. Besides, there exists a function ip, called wavelet 
such that {<0fc = 2 j / 2 ip(2 : > x — fc)} feeZ is a basis of Wj. 

In practise, for adaptive refinement we set the coarsest level jo and the 
finest level ji, jo < jx, and we decompose the space corresponding to the 
finest level on all the levels in between: 

V n = V j0 © W j0 © W i0+1 ®---®W h -x. 

A function / £ Vj 1 can then be decomposed as follows 

+00 ji — 1 +00 

where the (c^°); are the coefficients on the coarse mesh and the (dj)i the 
details at the different level in between. 



^1,^2+1 c 2fci+l,2fc 2 +2 ^1+1^2+1 



c j+1 




L 2fci+2,2A;2+1 



Fig. 1. Mesh refinement in 2D. 



In two dimensions, the prediction operator which defines the multireso- 
lution analysis is constructed by tensor product from the ID operator. In 
practise three different cases must be considered (see figure [1] for notations): 

1. Refinement in x (corresponding to points +1 2 k 2 anc ^ ^fci+i 2k 2 +2) : we 
use the ID prediction operator in x for fixed k 2 . 

2. Refinement in v (corresponding to points c^" 2fc 2 +i an< ^ c iti+2 2k 2 +i)' we 
use the ID prediction operator in v for fixed k\ . 

3. Refinement in v (corresponding to point 4t+i2): 2 +i) : we ^ rs ^ use the 
ID prediction operator in v for fixed k\ to determine the points which 
are necessary for applying the ID prediction operator in x for fixed 
which we then apply. 
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The corresponding wavelet bases are respectively of type i/)(x)ip(v), ip(x)ip(v) 
and ip(x)ip{v) where ip and ip are respectively the scaling function and the ID 
wavelet. We then obtain a 2D wavelet decomposition of the following form: 



3 The algorithms 

We want to numerically solve the Vlasov equation |T]) given an initial value 
of the distribution function Jo- 

We start by describing the method based on an interpolation using the 
wavelet decomposition of / in the non adaptive case. Then we overlay an 
adaptive algorithm to this method. 

For those two algorithms, we first pick the resolution levels for the phase- 
space meshes, from the coarsest jo to the hnest j\. Although these levels 
could be different in x and v, we consider here for the sake of conciseness and 
clarity that they are identical. 

We also compute our scaling function on a very fine grid so that we can 
obtain with enough precision its value at any point. 

3.1 The non adaptive algorithm 

We are working in this case on the finest level corresponding to ji keeping 
all the points. 

Initialization: We decompose the initial condition in the wavelet basis 
by computing the coefficients Ck lt k 2 of the decomposition in Vj for the coarse 
mesh, and then adding the details d? k k2 in the detail spaces Wj for all the 
other levels j = jo, . . . ,j\ — 1. We then compute the initial electric field. 

Time iterations: 

— Advection in x: We start by computing for each mesh point the origin 
of the corresponding characteristic exactly, the displacement being VjAt. 
As we do not necessarily land on a mesh point, we compute the values of 
the distribution function at the intermediate time level, denoted by /*, at 
the origin of the characteristics by interpolation from /". We use for this 
the wavelet decomposition applied to /" from which we can compute 
/" at any point in phase space. 




(3) 
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— Computation of the electric field: We compute the charge density 
by integrating /* with respect to v, then the electric field by solving the 
Poisson equation (this step vanishes for the linear case of the rotating 
cylinder where the advection field is exactly known). 

— Advection in v: We start by computing exactly the origin of the char- 
acteristic for each mesh point, the displacement being E(t n , Xi)At. As we 
do not necessarily land on a mesh point, we compute the values of the 
distribution function at the intermediate time level, denoted by f n+1 , at 
the origin of the characteristics by interpolation from /*. We use for this 
the wavelet decomposition of /* given by ([3]) used at the previous step. 

3.2 The adaptive algorithm 

In the initialization phase, we first compute the wavelet decomposition of 
the initial condition /o, and then proceed by compressing it, i.e. eliminating 
the details which are smaller than a threshold that we impose. We then 
construct an adaptive mesh which, from all the possible points at all the 
levels between our coarsest and finest, contains only those of the coarsest 
and those corresponding to details which are above the threshold. We denote 
by G this mesh. 

— Prediction in x: We predict the positions of points where the details 
should be important at the next time split step by advancing in x the 
characteristics originating from the points of the mesh G. For this we use an 
explicit Euler scheme for the numerical integration of the characteristics. 
Then we retain the grid points, at one level finer as the starting point, 
surrounding the end point the characteristic. 

— Construction of mesh G: From the predicted mesh G, we construct 
the mesh G where the values of the distribution at the next time step 
shall be computed. This mesh G contains exactly the points necessary for 
computing the wavelet transform of /* at the points of G. 

— Advection in x: As in the non adaptive case. 

— Wavelet transform of /*: We compute the Ck and dk coefficients at the 
points of G from the values of /* at the points of G. 

— Compression: We eliminate the points of G where the details dk are lower 
than the fixed threshold. 

— Computation of the electric field: As in the non adaptive case. 

— Prediction in v: As for the prediction in x. 

— Construction of mesh G: As previously. This mesh G contains exactly 
the points necessary for computing the wavelet transform of f n+1 at the 
points of G determined in the prediction in v step. 

— Advection in v. As in the non adaptive case. 

— Wavelet transform of /™ +1 : We compute the Ck and dk at the points of 
G from the values of f n+1 at the points of G. 

— Compression: We eliminate the points of G where the details dk are lower 
than the fixed threshold. 
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4 Numerical results 



We show here our first results obtained with the adaptive method. We con- 
sider first a linear problem, namely the test case of the rotating cylinder 
introduced by Zalesak [10] to test advection schemes. Then we consider a 
classical nonlinear Vlasov-Poisson test case, namely the two stream instabil- 
ity. 

4.1 The slit rotating cylinder 

We consider the following initial condition: 

1 if Vx 2 + v 2 < 0.5 and if x < or |v| > 0.125, 
else. 



f(0,x,v) 



The computational domain is [—0.5,0.5] x [—0.5,0.5]. 

The advection field is (v, —x), which corresponds to the Vlasov equation 
with an applied electric field E app (x, t) = —x and without self-consistent field. 
Figure [2] represents the evolution of the rotating cylinder on a half turn with 
a coarse mesh of 16 x 16 points and 4 adaptive refinement levels. We notice 
that the cylinder is well represented and that the mesh points concentrate 
along the discontinuities. 



4.2 The two-stream instability 

We consider two streams symmetric with respect to v = and represented 
by the initial distribution function 

/(0, x, v) = — v 2 exp(— v 2 /2)(1 + a cos(fco %)), 
\/2tt 

with a = 0.25, kg = 0.5, and L = 2ir/ko. We use a maximum of N x = 128 
points in the x direction, and N v = 128 points in the v direction with v ma x = 
7, and a time step At = 1/8. The solution varies first very slowly and then fine 
scales are generated. Between times of around t ~ 20 uo^ 1 and t ~ 40 w" 1 , 
the instability increases rapidly and a hole appears in the middle of the 
computational domain. After t = 45 lu^ 1 until the end of the simulation, 
particles inside the hole are trapped. On figure [3] we show a snapshot of the 
distribution function at times t = 5 lo^ 1 and t — 30 oj" 1 for a coarse mesh of 
16 x 16 points and 3 levels of refinement. The adaptive method reproduces 
well the results obtained in the non adaptive case. 



5 Conclusion 

In this paper we have described a new method for the numerical resolution of 
the Vlasov equation using an adaptive mesh of phase-space. The adaptive al- 
gorithm is based on a multiresolution analysis. It performs qualitatively well. 
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Fig. 2. Rotating cylinder: evolution for a coarse mesh of 2 4 x 2 4 points and 4 adap- 
tive refinement levels. Snapshots of the cylinder and the corresponding adaptive 
mesh: (upper) after one time step, (lower) after 1/2 turn. 




Fig. 3. Two stream instability for a coarse mesh of 2 4 x 2 4 , and 3 adaptive refine- 
ment levels, (left) at time t = hu^ 1 , (right) at time t = SOu)- 1 . 
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However, there is a large overhead due to the handling of the adaptive mesh 
which has not been optimized yet. The performance of the code needs to be 
improved before we can recommend this technique for actual computations. 
We are currently working on optimizing the code and trying different kinds 
of wavelets, as well as obtaining error estimates for the adaptive method. 
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